clc;
clear;

%% main主脚本
tic;
startTime = cputime;
pde = process_pde_t(1);
vector_h = [1/4 1/8 1/16];
iter_L = 5;
basis_type = "P2";

% process_showbasis(basis_type, "x");

u_Linf = zeros(size(vector_h,2),1);
u_L2 = zeros(size(vector_h,2),1);
u_H1 = zeros(size(vector_h,2),1);
p_Linf = zeros(size(vector_h,2),1);
p_L2 = zeros(size(vector_h,2),1);
p_H1 = zeros(size(vector_h,2),1);

for i = 1:size(vector_h,2)
    h = vector_h(i);
    mesh = process_mesh(pde, h, basis_type);

    % process_showmesh(mesh);

    O1 = sparse(mesh.Np, mesh.Np);
    O2 = sparse(mesh.Nb, mesh.Np);
    O3 = sparse(mesh.Nb, mesh.Nb);
    O4 = zeros(mesh.Np, 1);
    A1 = assemble_stiffness_matrix(mesh, mesh.basis_type, "dx", mesh.basis_type, "dx", pde.nu);
    A2 = assemble_stiffness_matrix(mesh, mesh.basis_type, "dy", mesh.basis_type, "dy", pde.nu);
    A3 = assemble_stiffness_matrix(mesh, mesh.basis_type, "dx", mesh.basis_type, "dy", pde.nu);
    A5 = assemble_stiffness_matrix(mesh, "P1", "x", mesh.basis_type, "dx", -1);
    A6 = assemble_stiffness_matrix(mesh, "P1", "x", mesh.basis_type, "dy", -1);
    A = [2*A1+A2, A3, A5;
        A3', 2*A2+A1, A6;
        A5', A6', O1];
    Me = assemble_mass_matrix(mesh);
    M = [Me, O3, O2;
        O3, Me, O2;
        O2', O2', O1];
    U1 = [pde.u1(mesh.node(:,1),mesh.node(:,2),0);
        pde.u2(mesh.node(:,1),mesh.node(:,2),0);
        pde.p(mesh.node(1:mesh.Np,1),mesh.node(1:mesh.Np,2),0)];
    dt = 8*h^3;
    Nt = (pde.domain(6) - pde.domain(5))/dt;
    for m = 0:Nt-1
        t2 = (m+1)*dt;
        [b1, b2] = assemble_load_vector(mesh, pde, t2);
        bt2 = [b1; b2; O4];
        U_old = U1;
        for l = 1:iter_L
            AN1 = assemble_nonlinear_matrix(mesh, "u1", "dx", mesh.basis_type, "x", mesh.basis_type, "x", U_old);
            AN2 = assemble_nonlinear_matrix(mesh, "u1", "x", mesh.basis_type, "dx", mesh.basis_type, "x", U_old);
            AN3 = assemble_nonlinear_matrix(mesh, "u2", "x", mesh.basis_type, "dy", mesh.basis_type, "x", U_old);
            AN4 = assemble_nonlinear_matrix(mesh, "u1", "dy", mesh.basis_type, "x", mesh.basis_type, "x", U_old);
            AN5 = assemble_nonlinear_matrix(mesh, "u2", "dx", mesh.basis_type, "x", mesh.basis_type, "x", U_old);
            AN6 = assemble_nonlinear_matrix(mesh, "u2", "dy", mesh.basis_type, "x", mesh.basis_type, "x", U_old);
            AN = [AN1+AN2+AN3, AN4, O2;
                AN5, AN6+AN2+AN3, O2;
                O2', O2', O1];
            bN1 = assemble_nonlinear_vector(mesh, "u1", "x", "u1", "dx", mesh.basis_type, "x", U_old);
            bN2 = assemble_nonlinear_vector(mesh, "u2", "x", "u1", "dy", mesh.basis_type, "x", U_old);
            bN3 = assemble_nonlinear_vector(mesh, "u1", "x", "u2", "dx", mesh.basis_type, "x", U_old);
            bN4 = assemble_nonlinear_vector(mesh, "u2", "x", "u2", "dy", mesh.basis_type, "x", U_old);
            bN = [bN1+bN2; bN3+bN4; O4];
            Al = M./dt + A + AN;
            bl = (M./dt)*U_old + bt2 + bN;
            [Al, bl] = process_boundary_t(Al, bl, mesh, pde, t2);
            Al(2*mesh.Nb+1,:) = 0;
            Al(2*mesh.Nb+1,2*mesh.Nb+1) = 1;
            bl(2*mesh.Nb+1,1) = pde.p(pde.domain(1), pde.domain(3), t2);
            U = Al\bl;
            U_old = U;
        end
        U1 = U;
    end

    process_plot(mesh, pde, U(1:mesh.Nb), basis_type);
    process_plot(mesh, pde, U(mesh.Nb+1:2*mesh.Nb), basis_type);
    process_plot(mesh, pde, U(2*mesh.Nb+1:end), "P1");
    process_field(mesh, U);

    [u_Linf(i,1), p_Linf(i,1)] = process_error_t(mesh, pde, U, "Linf");
    [u_L2(i,1), p_L2(i,1)] = process_error_t(mesh, pde, U, "L2");
    [u_H1(i,1), p_H1(i,1)] = process_error_t(mesh, pde, U, "H1");

    fprintf("已完成\t%d/%d\n", i, size(vector_h,2));
end

process_rate(u_Linf, vector_h);
process_rate(u_L2, vector_h);
process_rate(u_H1, vector_h);
process_rate(p_Linf, vector_h);
process_rate(p_L2, vector_h);
process_rate(p_H1, vector_h);

elapsedTOCTime = toc;
elapsedCPUTime = cputime - startTime;
disp(["TOC time(s)", num2str(elapsedTOCTime)]);
disp(["CPU time(s)", num2str(elapsedCPUTime)]);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%% 以下是脚本所用到的函数 %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% process_pde_t
function PDE = process_pde_t(index)
switch index
    case 1
        PDE.index = index;
        nu = 1;
        PDE.nu = nu;
        PDE.u1 = @(x,y,t) (x.^2.*y.^2 + exp(-y)).*cos(2*pi*t);
        PDE.u2 = @(x,y,t) (-(2/3)*x.*y.^3 + 2 - pi*sin(pi*x)).*cos(2*pi*t);
        PDE.p = @(x,y,t) -(2-pi*sin(pi*x)).*cos(2*pi*y).*cos(2*pi*t);
        PDE.u1dx = @(x,y,t) 2*x.*y.^2.*cos(2*pi*t);
        PDE.u1dy = @(x,y,t) -cos(2*pi*t).*(-2*y.*x.^2 + exp(-y));
        PDE.u2dx = @(x,y,t) -cos(2*pi*t).*((2*y.^3)./3 + pi^2*cos(x*pi));
        PDE.u2dy = @(x,y,t) -2*x.*y.^2.*cos(2*pi*t);
        PDE.pdx = @(x,y,t) pi^2*cos(2*pi*t).*cos(pi*x).*cos(2*pi*y);
        PDE.pdy = @(x,y,t) -2*pi*cos(2*pi*t).*sin(2*pi*y).*(pi*sin(pi*x) - 2);
        PDE.f1 = @(x,y,t) cos(2*pi*t).^2.*(-2*y.*x.^2 + exp(-y)).*((2*x.*y.^3)./3 + pi*sin(pi*x) - 2) - nu*(cos(2*pi*t).*(2*x.^2 + exp(-y)) + 2*y.^2.*cos(2*pi*t)) - 2*pi*sin(2*pi*t).*(exp(-y) + x.^2.*y.^2) + 2*x.*y.^2.*cos(2*pi*t).^2.*(exp(-y) + x.^2.*y.^2) + pi^2*cos(2*pi*t).*cos(pi*x).*cos(2*pi*y);
        PDE.f2 = @(x,y,t) 2*pi*sin(2*pi*t).*((2*x.*y.^3)./3 + pi*sin(pi*x) - 2) - cos(2*pi*t).^2.*((2*y.^3)./3 + pi^2*cos(x*pi)).*(exp(-y) + x.^2.*y.^2) - nu*(pi^3*cos(2*pi*t).*sin(pi*x) - 4*x.*y.*cos(2*pi*t)) - 2*pi*cos(2*pi*t).*sin(2*pi*y).*(pi*sin(pi*x) - 2) + 2*x.*y.^2.*cos(2*pi*t).^2.*((2*x.*y.^3)./3 + pi*sin(pi*x) - 2);
        PDE.domain = [0, 1, -0.25, 0, 0, 1];
        PDE.bdary = [0, 0, 0, 0];
    case 2
        PDE.index = index;
        nu = 1;
        PDE.nu = nu;
        PDE.u1 = @(x,y,t) (x.^2.*y.^2 + exp(-y)).*cos(2*pi*t);
        PDE.u2 = @(x,y,t) (-(2/3)*x.*y.^3 + 2 - pi*sin(pi*x)).*cos(2*pi*t);
        PDE.p = @(x,y,t) -(2-pi*sin(pi*x)).*cos(2*pi*y).*cos(2*pi*t);
        PDE.u1dx = @(x,y,t) 2*x.*y.^2.*cos(2*pi*t);
        PDE.u1dy = @(x,y,t) -cos(2*pi*t).*(-2*y.*x.^2 + exp(-y));
        PDE.u2dx = @(x,y,t) -cos(2*pi*t).*((2*y.^3)./3 + pi^2*cos(x*pi));
        PDE.u2dy = @(x,y,t) -2*x.*y.^2.*cos(2*pi*t);
        PDE.pdx = @(x,y,t) pi^2*cos(2*pi*t).*cos(pi*x).*cos(2*pi*y);
        PDE.pdy = @(x,y,t) -2*pi*cos(2*pi*t).*sin(2*pi*y).*(pi*sin(pi*x) - 2);
        PDE.f1 = @(x,y,t) cos(2*pi*t).^2.*(-2*y.*x.^2 + exp(-y)).*((2*x.*y.^3)./3 + pi*sin(pi*x) - 2) - nu*(cos(2*pi*t).*(2*x.^2 + exp(-y)) + 2*y.^2.*cos(2*pi*t)) - 2*pi*sin(2*pi*t).*(exp(-y) + x.^2.*y.^2) + 2*x.*y.^2.*cos(2*pi*t).^2.*(exp(-y) + x.^2.*y.^2) + pi^2*cos(2*pi*t).*cos(pi*x).*cos(2*pi*y);
        PDE.f2 = @(x,y,t) 2*pi*sin(2*pi*t).*((2*x.*y.^3)./3 + pi*sin(pi*x) - 2) - cos(2*pi*t).^2.*((2*y.^3)./3 + pi^2*cos(x*pi)).*(exp(-y) + x.^2.*y.^2) - nu*(pi^3*cos(2*pi*t).*sin(pi*x) - 4*x.*y.*cos(2*pi*t)) - 2*pi*cos(2*pi*t).*sin(2*pi*y).*(pi*sin(pi*x) - 2) + 2*x.*y.^2.*cos(2*pi*t).^2.*((2*x.*y.^3)./3 + pi*sin(pi*x) - 2);
        PDE.domain = [0, 1, -0.25, 0, 0, 1];
        PDE.bdary = [0, 1, 0, 0];
    case 3
        PDE.index = index;
        nu = 1;
        PDE.nu = nu;
        PDE.u1 = @(x,y,t) pi*sin(2*pi*y).*sin(pi*x).^2.*sin(t);
        PDE.u2 = @(x,y,t) -pi*sin(2*pi*x).*sin(pi*y).^2.*sin(t);
        PDE.p = @(x,y,t) 10*cos(pi*x).*sin(pi*y).*sin(t);
        PDE.u1dx = @(x,y,t) 2*pi^2*cos(pi*x).*sin(pi*x).*sin(2*pi*y).*sin(t);
        PDE.u1dy = @(x,y,t) 2*pi^2*cos(2*pi*y).*sin(pi*x).^2.*sin(t);
        PDE.u2dx = @(x,y,t) -2*pi^2*cos(2*pi*x).*sin(pi*y).^2.*sin(t);
        PDE.u2dy = @(x,y,t) -2*pi^2*cos(pi*y).*sin(2*pi*x).*sin(pi*y).*sin(t);
        PDE.pdx = @(x,y,t) -10*pi*sin(pi*x).*sin(pi*y).*sin(t);
        PDE.pdy = @(x,y,t) 10*pi*cos(pi*x).*cos(pi*y).*sin(t);
        PDE.f1 = @(x,y,t) pi*sin(pi*x).^2.*sin(2*pi*y).*cos(t) - nu*(2*pi^3*cos(pi*x).^2.*sin(2*pi*y).*sin(t) - 6*pi^3*sin(pi*x).^2.*sin(2*pi*y).*sin(t)) - 10*pi*sin(pi*x).*sin(pi*y).*sin(t) + 2*pi^3*cos(pi*x).*sin(pi*x).^3.*sin(2*pi*y).^2.*sin(t).^2 - 2*pi^3*cos(2*pi*y).*sin(pi*x).^2.*sin(2*pi*x).*sin(pi*y).^2.*sin(t).^2;
        PDE.f2 = @(x,y,t) nu*(2*pi^3*cos(pi*y).^2.*sin(2*pi*x).*sin(t) - 6*pi^3*sin(2*pi*x).*sin(pi*y).^2.*sin(t)) - pi*sin(2*pi*x).*sin(pi*y).^2.*cos(t) + 10*pi*cos(pi*x).*cos(pi*y).*sin(t) + 2*pi^3*cos(pi*y).*sin(2*pi*x).^2.*sin(pi*y).^3.*sin(t).^2 - 2*pi^3*cos(2*pi*x).*sin(pi*x).^2.*sin(pi*y).^2.*sin(2*pi*y).*sin(t).^2;
        PDE.domain = [0, 1, 0, 1, 0, 1];
        PDE.bdary = [0, 0, 0, 0];
    case 4
        PDE.index = index;
        nu = 1;
        PDE.nu = nu;
        PDE.u1 = @(x,y,t) pi*sin(2*pi*y).*sin(pi*x).^2.*sin(t);
        PDE.u2 = @(x,y,t) -pi*sin(2*pi*x).*sin(pi*y).^2.*sin(t);
        PDE.p = @(x,y,t) 10*cos(pi*x).*sin(pi*y).*sin(t);
        PDE.u1dx = @(x,y,t) 2*pi^2*cos(pi*x).*sin(pi*x).*sin(2*pi*y).*sin(t);
        PDE.u1dy = @(x,y,t) 2*pi^2*cos(2*pi*y).*sin(pi*x).^2.*sin(t);
        PDE.u2dx = @(x,y,t) -2*pi^2*cos(2*pi*x).*sin(pi*y).^2.*sin(t);
        PDE.u2dy = @(x,y,t) -2*pi^2*cos(pi*y).*sin(2*pi*x).*sin(pi*y).*sin(t);
        PDE.pdx = @(x,y,t) -10*pi*sin(pi*x).*sin(pi*y).*sin(t);
        PDE.pdy = @(x,y,t) 10*pi*cos(pi*x).*cos(pi*y).*sin(t);
        PDE.f1 = @(x,y,t) pi*sin(pi*x).^2.*sin(2*pi*y).*cos(t) - nu*(2*pi^3*cos(pi*x).^2.*sin(2*pi*y).*sin(t) - 6*pi^3*sin(pi*x).^2.*sin(2*pi*y).*sin(t)) - 10*pi*sin(pi*x).*sin(pi*y).*sin(t) + 2*pi^3*cos(pi*x).*sin(pi*x).^3.*sin(2*pi*y).^2.*sin(t).^2 - 2*pi^3*cos(2*pi*y).*sin(pi*x).^2.*sin(2*pi*x).*sin(pi*y).^2.*sin(t).^2;
        PDE.f2 = @(x,y,t) nu*(2*pi^3*cos(pi*y).^2.*sin(2*pi*x).*sin(t) - 6*pi^3*sin(2*pi*x).*sin(pi*y).^2.*sin(t)) - pi*sin(2*pi*x).*sin(pi*y).^2.*cos(t) + 10*pi*cos(pi*x).*cos(pi*y).*sin(t) + 2*pi^3*cos(pi*y).*sin(2*pi*x).^2.*sin(pi*y).^3.*sin(t).^2 - 2*pi^3*cos(2*pi*x).*sin(pi*x).^2.*sin(pi*y).^2.*sin(2*pi*y).*sin(t).^2;
        PDE.domain = [0, 1, 0, 1, 0, 1];
        PDE.bdary = [0, 1, 0, 0];
    case 5
        PDE.index = index;
        nu = 1;
        PDE.nu = nu;
        PDE.u1 = @(x,y,t) 10.*x.^2.*y.*cos(t).*(2.*y - 1).*(x - 1).^2.*(y - 1);
        PDE.u2 = @(x,y,t) -10.*x.*y.^2.*cos(t).*(2.*x - 1).*(x - 1).*(y - 1).^2;
        PDE.p = @(x,y,t) cos(t).*(20.*x - 10).*(2.*y - 1);
        PDE.u1dx = @(x,y,t) 20*x.*y.*cos(t).*(2*y - 1).*(x - 1).^2.*(y - 1) + 10*x.^2.*y.*cos(t).*(2*x - 2).*(2*y - 1).*(y - 1);
        PDE.u1dy = @(x,y,t) 10*x.^2.*cos(t).*(2*y - 1).*(x - 1).^2.*(y - 1) + 20*x.^2.*y.*cos(t).*(x - 1).^2.*(y - 1) + 10*x.^2.*y.*cos(t).*(2*y - 1).*(x - 1).^2;
        PDE.u2dx = @(x,y,t) -10*y.^2.*cos(t).*(2*x - 1).*(x - 1).*(y - 1).^2 - 20*x.*y.^2.*cos(t).*(x - 1).*(y - 1).^2 - 10*x.*y.^2.*cos(t).*(2*x - 1).*(y - 1).^2;
        PDE.u2dy = @(x,y,t) -20.*x.*y.*cos(t).*(2*x - 1).*(x - 1).*(y - 1).^2 - 10*x.*y.^2.*cos(t).*(2*x - 1).*(2*y - 2).*(x - 1);
        PDE.pdx = @(x,y,t) 20.*cos(t).*(2*y - 1);
        PDE.pdy = @(x,y,t) 2*cos(t).*(20*x - 10);
        PDE.f1 = @(x,y,t) 20*cos(t).*(2*y - 1) - nu*(40*x.^2.*cos(t).*(x - 1).^2.*(y - 1) + 20*x.^2.*cos(t).*(2*y - 1).*(x - 1).^2 + 40*x.^2.*y.*cos(t).*(x - 1).^2 + 20*x.^2.*y.*cos(t).*(2*y - 1).*(y - 1) + 20*y.*cos(t).*(2*y - 1).*(x - 1).^2.*(y - 1) + 40*x.*y.*cos(t).*(2*x - 2).*(2*y - 1).*(y - 1)) - 10*x.^2.*y.*sin(t).*(2*y - 1).*(x - 1).^2.*(y - 1) - 10*x.*y.^2.*cos(t).*(2*x - 1).*(x - 1).*(y - 1).^2.*(10*x.^2.*cos(t).*(2*y - 1).*(x - 1).^2.*(y - 1) + 20*x.^2.*y.*cos(t).*(x - 1).^2.*(y - 1) + 10*x.^2.*y.*cos(t).*(2.*y - 1).*(x - 1).^2) + 10*x.^2.*y.*cos(t).*(2*y - 1).*(20*x.*y.*cos(t).*(2*y - 1).*(x - 1).^2.*(y - 1) + 10*x.^2.*y.*cos(t).*(2*x - 2).*(2*y - 1).*(y - 1)).*(x - 1).^2.*(y - 1);
        PDE.f2 = @(x,y,t) 2*cos(t).*(20*x - 10) + nu*(40*y.^2.*cos(t).*(x - 1).*(y - 1).^2 + 20*y.^2.*cos(t).*(2*x - 1).*(y - 1).^2 + 40*x.*y.^2.*cos(t).*(y - 1).^2 + 20*x.*y.^2.*cos(t).*(2*x - 1).*(x - 1) + 20*x.*cos(t).*(2*x - 1).*(x - 1).*(y - 1).^2 + 40*x.*y.*cos(t).*(2*x - 1).*(2*y - 2).*(x - 1)) + 10*x.*y.^2.*sin(t).*(2*x - 1).*(x - 1).*(y - 1).^2 - 10*x.^2.*y.*cos(t).*(2*y - 1).*(x - 1).^2.*(y - 1).*(10*y.^2.*cos(t).*(2*x - 1).*(x - 1).*(y - 1).^2 + 20*x.*y.^2.*cos(t).*(x - 1).*(y - 1).^2 + 10*x.*y.^2.*cos(t).*(2*x - 1).*(y - 1).^2) + 10*x.*y.^2.*cos(t).*(2*x - 1).*(20*x.*y.*cos(t).*(2*x - 1).*(x - 1).*(y - 1).^2 + 10*x.*y.^2.*cos(t).*(2*x - 1).*(2*y - 2).*(x - 1)).*(x - 1).*(y - 1).^2;
        PDE.domain = [0, 1, 0, 1, 0, 1];
        PDE.bdary = [0, 0, 0, 0];
    case 6
        PDE.index = index;
        nu = 1;
        PDE.nu = nu;
        PDE.u1 = @(x,y,t) 10.*x.^2.*y.*cos(t).*(2.*y - 1).*(x - 1).^2.*(y - 1);
        PDE.u2 = @(x,y,t) -10.*x.*y.^2.*cos(t).*(2.*x - 1).*(x - 1).*(y - 1).^2;
        PDE.p = @(x,y,t) cos(t).*(20.*x - 10).*(2.*y - 1);
        PDE.u1dx = @(x,y,t) 20*x.*y.*cos(t).*(2*y - 1).*(x - 1).^2.*(y - 1) + 10*x.^2.*y.*cos(t).*(2*x - 2).*(2*y - 1).*(y - 1);
        PDE.u1dy = @(x,y,t) 10*x.^2.*cos(t).*(2*y - 1).*(x - 1).^2.*(y - 1) + 20*x.^2.*y.*cos(t).*(x - 1).^2.*(y - 1) + 10*x.^2.*y.*cos(t).*(2*y - 1).*(x - 1).^2;
        PDE.u2dx = @(x,y,t) -10*y.^2.*cos(t).*(2*x - 1).*(x - 1).*(y - 1).^2 - 20*x.*y.^2.*cos(t).*(x - 1).*(y - 1).^2 - 10*x.*y.^2.*cos(t).*(2*x - 1).*(y - 1).^2;
        PDE.u2dy = @(x,y,t) -20.*x.*y.*cos(t).*(2*x - 1).*(x - 1).*(y - 1).^2 - 10*x.*y.^2.*cos(t).*(2*x - 1).*(2*y - 2).*(x - 1);
        PDE.pdx = @(x,y,t) 20.*cos(t).*(2*y - 1);
        PDE.pdy = @(x,y,t) 2*cos(t).*(20*x - 10);
        PDE.f1 = @(x,y,t) 20*cos(t).*(2*y - 1) - nu*(40*x.^2.*cos(t).*(x - 1).^2.*(y - 1) + 20*x.^2.*cos(t).*(2*y - 1).*(x - 1).^2 + 40*x.^2.*y.*cos(t).*(x - 1).^2 + 20*x.^2.*y.*cos(t).*(2*y - 1).*(y - 1) + 20*y.*cos(t).*(2*y - 1).*(x - 1).^2.*(y - 1) + 40*x.*y.*cos(t).*(2*x - 2).*(2*y - 1).*(y - 1)) - 10*x.^2.*y.*sin(t).*(2*y - 1).*(x - 1).^2.*(y - 1) - 10*x.*y.^2.*cos(t).*(2*x - 1).*(x - 1).*(y - 1).^2.*(10*x.^2.*cos(t).*(2*y - 1).*(x - 1).^2.*(y - 1) + 20*x.^2.*y.*cos(t).*(x - 1).^2.*(y - 1) + 10*x.^2.*y.*cos(t).*(2.*y - 1).*(x - 1).^2) + 10*x.^2.*y.*cos(t).*(2*y - 1).*(20*x.*y.*cos(t).*(2*y - 1).*(x - 1).^2.*(y - 1) + 10*x.^2.*y.*cos(t).*(2*x - 2).*(2*y - 1).*(y - 1)).*(x - 1).^2.*(y - 1);
        PDE.f2 = @(x,y,t) 2*cos(t).*(20*x - 10) + nu*(40*y.^2.*cos(t).*(x - 1).*(y - 1).^2 + 20*y.^2.*cos(t).*(2*x - 1).*(y - 1).^2 + 40*x.*y.^2.*cos(t).*(y - 1).^2 + 20*x.*y.^2.*cos(t).*(2*x - 1).*(x - 1) + 20*x.*cos(t).*(2*x - 1).*(x - 1).*(y - 1).^2 + 40*x.*y.*cos(t).*(2*x - 1).*(2*y - 2).*(x - 1)) + 10*x.*y.^2.*sin(t).*(2*x - 1).*(x - 1).*(y - 1).^2 - 10*x.^2.*y.*cos(t).*(2*y - 1).*(x - 1).^2.*(y - 1).*(10*y.^2.*cos(t).*(2*x - 1).*(x - 1).*(y - 1).^2 + 20*x.*y.^2.*cos(t).*(x - 1).*(y - 1).^2 + 10*x.*y.^2.*cos(t).*(2*x - 1).*(y - 1).^2) + 10*x.*y.^2.*cos(t).*(2*x - 1).*(20*x.*y.*cos(t).*(2*x - 1).*(x - 1).*(y - 1).^2 + 10*x.*y.^2.*cos(t).*(2*x - 1).*(2*y - 2).*(x - 1)).*(x - 1).*(y - 1).^2;
        PDE.domain = [0, 1, 0, 1, 0, 1];
        PDE.bdary = [0, 1, 0, 0];
    otherwise
        error("Invalid PDE index.");
end
fprintf("PDE index: %d\n", PDE.index);
fprintf("PDE domain: xmin,xmax,ymin,ymax,tmix,tmax\n");
disp(PDE.domain);
fprintf("PDE boundary: bottom,right,top,left (0=Dirichlet,1=Neumann,2=Robin)\n");
disp(PDE.bdary);
end

%% assemble_stiffness_matrix
function matrix_A = assemble_stiffness_matrix(mesh, B1_type, D1_type, B2_type, D2_type, cor)
basis_map = containers.Map(...
    {'P1', 'P1b', 'P2'},...
    {struct('Nb', mesh.Np, 'Nlb',3),...
    struct('Nb', mesh.Nb, 'Nlb',4),...
    struct('Nb', mesh.Nb, 'Nlb',6)});
B1 = basis_map(B1_type);
Nb1 = B1.Nb();
Nlb1 = B1.Nlb();
B2 = basis_map(B2_type);
Nb2 = B2.Nb();
Nlb2 = B2.Nlb();
rows = zeros(2*Nlb1*Nlb2,1);
cols = rows;
values = rows;
flag = 1;
for n = 1:mesh.Ne
    E = mesh.node(mesh.elem(n,:),:);
    V = E(1:3,:);
    [gauss_points_loc, gauss_weights_loc] = gauss_lengdre_loc(V, "2d");
    phi1 = basis_function_loc(V, gauss_points_loc(:,1), gauss_points_loc(:,2), B1_type, D1_type);
    phi2 = basis_function_loc(V, gauss_points_loc(:,1), gauss_points_loc(:,2), B2_type, D2_type);
    for i = 1:Nlb2
        for j = 1:Nlb1
            f = cor * phi1(:,j) .* phi2(:,i);
            r = gauss_weights_loc * f;
            rows(flag) = mesh.elem(n,i);
            cols(flag) = mesh.elem(n,j);
            values(flag) = r;
            flag = flag + 1;
        end
    end
end
matrix_A = sparse(rows, cols, values, Nb2, Nb1);
end

%% assemble_mass_matrix
function matrix_M = assemble_mass_matrix(mesh)
Nb1 = mesh.Nb;
Nlb1 = mesh.Nlb;
Nb2 = Nb1;
Nlb2 = Nlb1;
rows = zeros(2*Nlb1*Nlb2,1);
cols = rows;
values = rows;
flag = 1;
for n = 1:mesh.Ne
    E = mesh.node(mesh.elem(n,:),:);
    V = E(1:3,:);
    [gauss_points_loc, gauss_weights_loc] = gauss_lengdre_loc(V, "2d");
    phi = basis_function_loc(V, gauss_points_loc(:,1), gauss_points_loc(:,2), mesh.basis_type, "x");
    for i = 1:Nlb2
        for j = 1:Nlb1
            f = phi(:,j) .* phi(:,i);
            r = gauss_weights_loc * f;
            rows(flag) = mesh.elem(n,i);
            cols(flag) = mesh.elem(n,j);
            values(flag) = r;
            flag = flag + 1;
        end
    end
end
matrix_M = sparse(rows, cols, values, Nb2, Nb1);
end

%% assemble_load_vector
function [vector_b1, vector_b2] = assemble_load_vector(mesh, pde, t)
vector_b1 = zeros(mesh.Nb,1);
vector_b2 = vector_b1;
for n = 1:mesh.Ne
    E = mesh.node(mesh.elem(n,:),:);
    V = E(1:3,:);
    [gauss_points_loc, gauss_weights_loc] = gauss_lengdre_loc(V, "2d");
    phi = basis_function_loc(V, gauss_points_loc(:,1), gauss_points_loc(:,2), mesh.basis_type, "x");
    ff1 = pde.f1(gauss_points_loc(:,1), gauss_points_loc(:,2), t);
    ff2 = pde.f2(gauss_points_loc(:,1), gauss_points_loc(:,2), t);
    for i = 1:mesh.Nlb
        f1 = ff1 .* phi(:,i);
        f2 = ff2 .* phi(:,i);
        r1 = gauss_weights_loc * f1;
        r2 = gauss_weights_loc * f2;
        vector_b1(mesh.elem(n,i),1) = vector_b1(mesh.elem(n,i),1) + r1;
        vector_b2(mesh.elem(n,i),1) = vector_b2(mesh.elem(n,i),1) + r2;
    end
end
end

%% assemble_nonlinear_matrix
function matrix_A = assemble_nonlinear_matrix(mesh, u_type, du_type, B1_type, D1_type, B2_type, D2_type, U_old)
basis_map = containers.Map(...
    {'P1', 'P1b', 'P2'},...
    {struct('Nb', mesh.Np, 'Nlb',3),...
    struct('Nb', mesh.Nb, 'Nlb',4),...
    struct('Nb', mesh.Nb, 'Nlb',6)});
B1 = basis_map(B1_type);
Nb1 = B1.Nb();
Nlb1 = B1.Nlb();
B2 = basis_map(B2_type);
Nb2 = B2.Nb();
Nlb2 = B2.Nlb();
u_map = containers.Map(...
    {'u1', 'u2'},...
    {struct('shift', 0),...
    struct('shift', mesh.Nb)});
U = u_map(u_type);
shift = U.shift;
rows = zeros(2*Nlb1*Nlb2,1);
cols = rows;
values = rows;
flag = 1;
for n = 1:mesh.Ne
    E = mesh.node(mesh.elem(n,:),:);
    V = E(1:3,:);
    [gauss_points_loc, gauss_weights_loc] = gauss_lengdre_loc(V, "2d");
    phiu = basis_function_loc(V, gauss_points_loc(:,1), gauss_points_loc(:,2), "P1", du_type);
    phi1 = basis_function_loc(V, gauss_points_loc(:,1), gauss_points_loc(:,2), B1_type, D1_type);
    phi2 = basis_function_loc(V, gauss_points_loc(:,1), gauss_points_loc(:,2), B2_type, D2_type);
    item = phiu * U_old(mesh.elem(n,1:3)+shift);
    for i = 1:Nlb2
        for j = 1:Nlb1
            f = item .* phi1(:,j) .* phi2(:,i);
            r = gauss_weights_loc * f;
            rows(flag) = mesh.elem(n,i);
            cols(flag) = mesh.elem(n,j);
            values(flag) = r;
            flag = flag + 1;
        end
    end
end
matrix_A = sparse(rows, cols, values, Nb2, Nb1);
end

%% assemble_nonlinear_vector
function vector_b = assemble_nonlinear_vector(mesh, u1_type, du1_type, u2_type, du2_type, B_type, D_type, U_old)
u_map = containers.Map(...
    {'u1', 'u2'},...
    {struct('shift', 0),...
    struct('shift', mesh.Nb)});
U1 = u_map(u1_type);
shift1 = U1.shift;
U2 = u_map(u2_type);
shift2 = U2.shift;
vector_b = zeros(mesh.Nb,1);
for n = 1:mesh.Ne
    E = mesh.node(mesh.elem(n,:),:);
    V = E(1:3,:);
    [gauss_points_loc, gauss_weights_loc] = gauss_lengdre_loc(V, "2d");
    phi1 = basis_function_loc(V, gauss_points_loc(:,1), gauss_points_loc(:,2), "P1", du1_type);
    phi2 = basis_function_loc(V, gauss_points_loc(:,1), gauss_points_loc(:,2), "P1", du2_type);
    phi3 = basis_function_loc(V, gauss_points_loc(:,1), gauss_points_loc(:,2), B_type, D_type);
    item1 = phi1 * U_old(mesh.elem(n,1:3)+shift1);
    item2 = phi2 * U_old(mesh.elem(n,1:3)+shift2);
    for j = 1:mesh.Nlb
        f =  item1 .* item2 .* phi3(:,j);
        r = gauss_weights_loc * f;
        vector_b(mesh.elem(n,j),1) = vector_b(mesh.elem(n,j),1) + r;
    end
end
end

%% process_boundary
function [A, b] = process_boundary_t(A, b, mesh, pde, t)
shift_u2 = mesh.Nb;
% Neumann
vector_b = zeros(2*mesh.Nb+mesh.Np,1);
for k = 1:size(pde.bdary,2)
    if 1 == pde.bdary(1,k)
        bde = boundary_edges(mesh, k);
        for n = 1:size(bde,1)
            points = mesh.node(bde(n,2:3),:);
            AB = [points(2,1)-points(1,1); points(2,2)-points(1,2)];
            vector_n = [0, 1; -1, 0] * AB;
            vector_n = vector_n./norm(vector_n);
            fp1 = @(x,y,t) (2*pde.nu*pde.u1dx(x,y,t)-pde.p(x,y,t))*vector_n(1) + pde.nu*(pde.u1dy(x,y,t)+pde.u2dx(x,y,t))*vector_n(2);
            fp2 = @(x,y,t) pde.nu*(pde.u1dy(x,y,t)+pde.u2dx(x,y,t))*vector_n(1) + (2*pde.nu*pde.u2dy(x,y,t)-pde.p(x,y,t))*vector_n(2);
            nk = bde(n,1);
            E = mesh.node(mesh.elem(nk,:),:);
            V = E(1:3,:);
            [gauss_points_loc, gauss_weights_loc] = gauss_lengdre_loc(points, "1d");
            phi = basis_function_loc(V, gauss_points_loc(:,1), gauss_points_loc(:,2), mesh.basis_type, "x");
            ffp1 = fp1(gauss_points_loc(:,1), gauss_points_loc(:,2), t);
            ffp2 = fp2(gauss_points_loc(:,1), gauss_points_loc(:,2), t);
            for i = 1:mesh.Nlb
                f1 = ffp1 .* phi(:,i);
                f2 = ffp2 .* phi(:,i);
                r1 = gauss_weights_loc * f1;
                r2 = gauss_weights_loc * f2;
                vector_b(mesh.elem(nk,i),1) = vector_b(mesh.elem(nk,i),1) + r1;
                vector_b(mesh.elem(nk,i)+shift_u2,1) = vector_b(mesh.elem(nk,i)+shift_u2,1) + r2;
            end
        end
    end
end
b = b + vector_b;
% Dirichlet
for k = 1:size(pde.bdary,2)
    if 0 == pde.bdary(1,k)
        bdn1 = boundary_nodes(mesh, k);
        A(bdn1,:) = 0;
        A(sub2ind(size(A), bdn1, bdn1)) = 1;
        b(bdn1) = pde.u1(mesh.node(bdn1,1), mesh.node(bdn1,2), t);
        bdn2 = bdn1 + shift_u2;
        A(bdn2,:) = 0;
        A(sub2ind(size(A), bdn2, bdn2)) = 1;
        b(bdn2) = pde.u2(mesh.node(bdn1,1), mesh.node(bdn1,2), t);
    end
end
end

%% process_error
function [error_u, error_p] = process_error_t(mesh, pde, U, error_type)
t_end = pde.domain(6);
shift_u2 = mesh.Nb;
shift_p = 2*mesh.Nb;
switch error_type
    case "Linf"
        err_u1 = 0;
        err_u2 = 0;
        err_p = 0;
        for n = 1:mesh.Ne
            E = mesh.node(mesh.elem(n,:),:);
            V = E(1:3,:);
            [gauss_points_loc, ~] = gauss_lengdre_loc(V, "2d");
            phi = basis_function_loc(V, gauss_points_loc(:,1), gauss_points_loc(:,2), mesh.basis_type, "x");
            psi = basis_function_loc(V, gauss_points_loc(:,1), gauss_points_loc(:,2), "P1", "x");
            u1h = phi * U(mesh.elem(n,:));
            u2h = phi * U(mesh.elem(n,:)+shift_u2);
            ph = psi * U(mesh.elem(n,1:3)+shift_p);
            f1 = max(abs(pde.u1(gauss_points_loc(:,1),gauss_points_loc(:,2),t_end) - u1h));
            f2 = max(abs(pde.u2(gauss_points_loc(:,1),gauss_points_loc(:,2),t_end) - u2h));
            f3 = max(abs(pde.p(gauss_points_loc(:,1),gauss_points_loc(:,2),t_end) - ph));
            if f1 > err_u1
                err_u1 = f1;
            end
            if f2 > err_u2
                err_u2 = f2;
            end
            if f3 > err_p
                err_p = f3;
            end
        end
        error_u = max(err_u1, err_u2);
        error_p = err_p;
    case "L2"
        err_u1 = 0;
        err_u2 = 0;
        err_p = 0;
        for n = 1:mesh.Ne
            E = mesh.node(mesh.elem(n,:),:);
            V = E(1:3,:);
            [gauss_points_loc, gauss_weights_loc] = gauss_lengdre_loc(V, "2d");
            phi = basis_function_loc(V, gauss_points_loc(:,1), gauss_points_loc(:,2), mesh.basis_type, "x");
            psi = basis_function_loc(V, gauss_points_loc(:,1), gauss_points_loc(:,2), "P1", "x");
            u1h = phi * U(mesh.elem(n,:));
            u2h = phi * U(mesh.elem(n,:)+shift_u2);
            ph = psi * U(mesh.elem(n,1:3)+shift_p);
            f1 = (pde.u1(gauss_points_loc(:,1),gauss_points_loc(:,2),t_end) - u1h).^2;
            f2 = (pde.u2(gauss_points_loc(:,1),gauss_points_loc(:,2),t_end) - u2h).^2;
            f3 = (pde.p(gauss_points_loc(:,1),gauss_points_loc(:,2),t_end) - ph).^2;
            r1 = gauss_weights_loc * f1;
            r2 = gauss_weights_loc * f2;
            r3 = gauss_weights_loc * f3;
            err_u1 = err_u1 + r1;
            err_u2 = err_u2 + r2;
            err_p = err_p + r3;
        end
        error_u = sqrt(err_u1 + err_u2);
        error_p = sqrt(err_p);
    case "H1"
        err_u1x = 0;
        err_u1y = 0;
        err_u2x = 0;
        err_u2y = 0;
        err_px = 0;
        err_py = 0;
        for n = 1:mesh.Ne
            E = mesh.node(mesh.elem(n,:),:);
            V = E(1:3,:);
            [gauss_points_loc, gauss_weights_loc] = gauss_lengdre_loc(V, "2d");
            phix = basis_function_loc(V, gauss_points_loc(:,1), gauss_points_loc(:,2), mesh.basis_type, "dx");
            phiy = basis_function_loc(V, gauss_points_loc(:,1), gauss_points_loc(:,2), mesh.basis_type, "dy");
            psix = basis_function_loc(V, gauss_points_loc(:,1), gauss_points_loc(:,2), "P1", "dx");
            psiy = basis_function_loc(V, gauss_points_loc(:,1), gauss_points_loc(:,2), "P1", "dy");
            u1hx = phix * U(mesh.elem(n,:));
            u1hy = phiy * U(mesh.elem(n,:));
            u2hx = phix * U(mesh.elem(n,:)+shift_u2);
            u2hy = phiy * U(mesh.elem(n,:)+shift_u2);
            phx = psix * U(mesh.elem(n,1:3)+shift_p);
            phy = psiy * U(mesh.elem(n,1:3)+shift_p);
            f1x = (pde.u1dx(gauss_points_loc(:,1),gauss_points_loc(:,2),t_end) - u1hx).^2;
            f1y = (pde.u1dy(gauss_points_loc(:,1),gauss_points_loc(:,2),t_end) - u1hy).^2;
            f2x = (pde.u2dx(gauss_points_loc(:,1),gauss_points_loc(:,2),t_end) - u2hx).^2;
            f2y = (pde.u2dy(gauss_points_loc(:,1),gauss_points_loc(:,2),t_end) - u2hy).^2;
            f3x = (pde.pdx(gauss_points_loc(:,1),gauss_points_loc(:,2),t_end) - phx).^2;
            f3y = (pde.pdy(gauss_points_loc(:,1),gauss_points_loc(:,2),t_end) - phy).^2;
            r1x = gauss_weights_loc * f1x;
            r1y = gauss_weights_loc * f1y;
            r2x = gauss_weights_loc * f2x;
            r2y = gauss_weights_loc * f2y;
            r3x = gauss_weights_loc * f3x;
            r3y = gauss_weights_loc * f3y;
            err_u1x = err_u1x + r1x;
            err_u1y = err_u2y + r1y;
            err_u2x = err_u2x + r2x;
            err_u2y = err_u2y + r2y;
            err_px = err_px + r3x;
            err_py = err_py + r3y;
        end
        error_u = sqrt(err_u1x + err_u1y + err_u2x + err_u2y);
        error_p = sqrt(err_px + err_py);
    otherwise
        error("Invalid error type.");
end
end